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Abstract 

In this paper, we characterize the network stabihty region (capacity region) of multi-queue multi- 
server (MQMS) queueing systems with stationary channel distribution and stationary arrival processes. 
The stability region is specified by a finite set of linear inequalities. We first show that the stability 
region is a poly tope characterized by the finite set of its facet defining hyperplanes. We explicitly 
determine the coefficients of the linear inequalities describing the facet defining hyperplanes of the 
stability region polytope. We further derive the necessary and sufficient conditions for the stability of 
the system for general arrival processes with finite first and second moments. For the case of stationary 
arrival processes, the derived conditions characterize the system stability region. Furthermore, we obtain 
an upper bound for the average queueing delay of Maximum Weight (MW) server allocation policy 
which has been shown in the literature to be a throughput optimal poHcy for MQMS systems. Using a 
similar approach, we can characterize the stability region for a fluid model MQMS system. However, 
the stability region of the fluid model system is described by an infinite number of linear inequalities 
since in this case the stability region is a convex surface. We present an example where we show that in 
some cases depending on the channel distribution, the stability region can be characterized by a finite 
set of non-linear inequalities instead of an infinite number of linear inequalities. 

I. Introduction 

Optimal stochastic network control is one of the primary goals in the design of emerging 
wireless networks. One of the objectives of stochastic control in wireless networks is to enable 
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cross layer designs to achieve stochastically optimal resource allocation in the physical and MAC 
layers, coupled with flow control/utility optimization strategies in transport layer and routing in 
the network layer. Examples of network resources at the MAC layer of wireless systems are 
OFDM subcarriers and CDMA codes and at the physical layer is transmission power. Apart 
from the resource allocation problem in physical and MAC layers, flow control plays a crucial 
role in ensuring the system stability while achieving a level of network fairness. A flow control 
strategy must decide how much fraction of the injected traffic to the system must be admitted to 
assure the stability of queues in the network layer while achieving the optimal network fairness 
among users. A quantitative measure of fairness which is widely used in literature (e.g. [[l|-[[9j) 
is to define a set of utility functions which illustrate the grade of satisfaction for each user 
n while it transmits data traffic at rate r. Consider a general network (wired or wireless) with 
source nodes. Suppose that A^, n = 1, 2, be the traffic arrival rate of each user n. For such 
a network, the stability region or network capacity region is defined as the closure of the set of 
all arrival rate vectors for which there exists a resource allocation policy that can stabilize the 
system [|T|, Q, pO| . Let us denote the network stability region by A. The flow control strategy 



determines the admitted rates from each user n by solving the following flow maximization 
problem. 

N 

Maximize: ^fn{rn) (1) 

n=l 

Subject to: ^ = (^i, ^2, ^at) ^ A 

< < Vn = 1,2, A^ 

By solving the above maximization problem, each node adjusts its admitted traffic rate to the 
computed optimal point from ([!]). We assume that functions fni^n) are non-decreasing and 
concave. The choice of functions fni'^n) depends on the desired fairness properties in the network. 
For example, choosing the utility function /(r) = log(l + ^r) (for some large constant parameter 
;5 > 0) will result in proportional fairness behaviour Q, Q. 

To solve the optimization problem in ([T]), we need to have the stability region of the system. 
Note that the stability region is unique for each network and independent of the resource 
allocation policy. Although stability region has been described and characterized for a general 
network in [[T|, Q, [[10), JTT| as given by the convex hull of a set of fixed points, such an 
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implicit characterization cannot be directly used to define the constraints of network optimization 
problem ([!]). The main goal of this paper is to introduce a linear algebraic characterization of 
the stability region for stationary Multi-Queue Multi-Server (MQMS) queueing systems. Such 
queueing systems can be used to model practical multi-user wireless networks with multiple 
orthogonal sub-channels such as OFDM sub-carriers. In such networks, there is a set of users 
generating random packet arrivals and a set of shared orthogonal sub-channels (servers) that are 
assigned to users according to certain rules. Because of users mobility, environmental changes, 
fading, etc., the channel quality of each user to each server is changing randomly with time. 
Therefore, resource allocation in such networks can be modeled as a server allocation problem in 
multi-queue multi-server queueing systems with time varying channel conditions [[T2)-[[T8|. Our 
focus in this paper would be on the server allocation problem in MQMS systems with stationary 
channel distribution for which we will introduce a linear algebraic representation of the network 
stability region (i.e.. A). Specifically, we will determine explicitly all the coefficients of the linear 
inequalities that describe the stability region. These inequalities then can be tabulated and used 
as the constraints of flow optimization problems similar to ([T]). 

The stability problem in wireless queueing networks was mainly addressed in [[T|, Q, pO| , 



[ [TT| , [[19). In [ [T0| |, authors introduced the notion of stability region of a queueing network. 



They considered a time slotted system in their work and assumed that arrival processes are 
i.i.d. sequences and the queue length process is a Markov process. In [^], they characterized 
the network stability region of multi-queue single-server systems with time varying ON-OFF 
connectivities. They also proved that for a symmetric system (with the same arrival and con- 
nectivity statistics for all the queues), LCQ (Longest Connected Queue) policy maximizes the 
stability region and also provides the optimal performance in terms of average queue occupancy 
(or equivalently average queueing delay). In [[T|, Q and |[20|, the notion of network stability 



region of a wireless network was introduced for more general arrival and queue length processes. 
Furthermore, Lyapunov drift techniques were applied in [|T) and Q to analyze the stability of 
the proposed policies for stochastic optimization problems in wireless networks. 

The problem of server allocation in multi-queue multi- server systems with time varying 
connectivities was mainly addressed in [[T2|, [[13]], [[TSj, [[T6|, [[T8|. In [ [T3] |, Maximum Weight 



(MW) policy was proposed as a throughput optimal server allocation policy for MQMS systems 
with stationary channel process. However, in [ |13[ the conditions on the arrival traffic to guarantee 
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the stability of MW were not explicitly mentioned. In our previous work in p8| |, we characterized 
the network stability region of multi-queue multi-server systems with time varying ON-OFF 
channels. We also obtained an upper bound for the average queueing delay of AS/LCQ (Any 
Server/Longest Connected Queue) policy which is the throughput optimal server allocation policy 
for such systems. 



References [ |T2] |, [[T4)-|[T6| study the optimal server allocation problem in terms of average 



queueing delay. In p2] , p4] , p31 l, the authors argue that in general, achieving instantaneous 
throughput and load balancing is impossible in a general MQMS system. However, they showed 
that this goal is attainable in the special case with ON-OFF channel processes. They also 
introduced MTLB (Maximum-Throughput Load-Balancing) policy and showed that this policy is 
minimizing a class of cost functions including total average delay for the case of two symmetric 
queues. The work in [ |16| | considers this problem for general number of symmetric queues and 



servers. Authors in p6| | characterized a class of Most Balancing (MB) policies among all work 
conserving policies which are minimizing, in stochastic ordering sense, a class of cost functions 
including total average delay. They used stochastic ordering and dynamic coupling arguments 
to show the optimality of MB policies for symmetric systems. 

In this paper, we will characterize the stability region of multi-queue multi- server queueing 
system with stationary channel and arrival processes. Toward this, the necessary and sufficient 
conditions for the stability of the system are derived for general arrival processes with finite first 
and second moments. For stationary arrival processes, these conditions establish the network 
stability region of such systems. Our contribution in this work is to characterize the stability 
region as a convex polytope specified by a finite set of linear inequalities that can be numerically 
tabulated and used to solve network optimization problems similar to ([T]) (Refer to Lemmas [3] to 
[7] and Theorems [T] and [2]). Later in the paper, we further introduce an upper bound for the average 



queueing delay of MW policy p3] | which is a throughput optimal policy for MQMS systems. 
We also study the stability of fluid model MQMS systems for which using the same approach 
as what we use in the packetized system, we characterize the linear algebraic representation of 
the stability region. We show that the stability region in this case is characterized by an infinite 
number of linear inequalities. However, by an example we show that depending on the channel 
distribution and the dimension of the system, we may characterize the stability region by a 
limited number of non-linear inequalities instead of infinite number of linear inequalities. 
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The rest of this paper is organized as follows. In section [11} we introduce the notation required 
throughout this paper. Section [Till describes the queueing model we focused on. In section [IV} we 



discuss about the notion and definition of strong stability in queueing networks. We also briefly 
elaborate on the Lyapunov drift technique used to prove the system stability [[T|, [[2). Moreover, 
we review some fundamental properties of polytopes in section |IV} In section |V} we will derive 



necessary and sufficient conditions for the stability of our model. We determine the coefficients 
of the linear inequalities describing the facets of the stability region polytope. We also find an 
upper bound for the average queue occupancy (or average queueing delay). Finally, we introduce 
fluid model MQMS systems and study the stability region of such queueing systems. Section 



VI presents the conclusions of our work. 



II. Notation 

In this section, we introduce basic notation used throughout the paper. Additional notation 
will be introduced when necessary. All the vectors are considered to be row vectors. By 1^ 
(0^), we denote a row vector of size K whose elements are all identically equal to "1" ("0"). 
The time average of a function f{t) is denoted by /(t), i.e., f{t) = lim^^oo ^ ^r=i fi^)- The 
operator " ® " is used for entry-wise multiplication of two matrices. The expectation of random 
processes (or random variables) is denoted by E[-]. The cardinality of a set is denoted by | • |. 
The operator for inner product of two vectors is (•,•). The boundary of a set is represented by 
hound{'). For any vector a = (o^i, 0^2, <^Ar) and a non-empty ordered (sub)set of indices U = 
{ui,U2, ^ {1,2, A^} and ui < U2 < ... < u\u\, we define au = {au^.au^, 

III. Model Description 

We consider a time slotted queueing system with equal length time slots and equal length 
packets. The model consists of a set of parallel queues A/" = {1, 2, A^} and a set of identical 
servers /C = {1, 2, K}. Each server can serve at most one queue at each time slot, i.e., we do 
not allow server sharing in the system. At each time slot t, the capacity of the link between each 
queue n G A/" and server A; G /C is assumed to be Cn^k{t) packets/time slot (Figure [T]), i.e., during 
time slot t, server k can serve at most Cn^k{t) packets of queue n successfully if it is allocated to 
queue n at that time slot. We assume that Cn^k{t) ^ where M = {m G Z+ | m < M}, for a 
given M. Therefore, at each time slot t, the channel state may be expressed by an x X matrix 
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Fig. 1: Multi-queue multi-server queueing system with stationary channel distribution 



C{t) = {Cn,k{t)) ,n G A/", G ]C,Cn,k{t) G M. The channel process is defined as {C'(t)}^^ 
with the state space S. Note that 5 is a finite set with \S\ = (M + 1)^^. We will label each 
element of 5 by a positive integer index s G {1,2,..., (M + 1)^^}. Suppose that the channel 
state matrix associated to the channel state s is denoted by Cg. The channel process is assumed 
to have stationary distribution with stationary probabilities tt^ = Pr(C(t) = Cg). 

Queues are fed by exogenous packet arrival processes yl^(t), n = 1, 2, A^, i.e., the number 
of packet arrivals to queue n during time slot t is represented by An{t). The arrival vector at 
time slot t is denoted by A{t) = {Ai{t), A2{t), AN{t)). For these processes, suppose that 
E[Al{t)] < A'^^^ < oc for all t. We assume that each queue has an infinite buffer space. 
We also assume that new arrivals are added to each queue at the end of each time slot. Let 
X{t) = {Xi{t), X^it)) be the queue length vector at the end of time slot t after adding new 
arrivals to the queues ^ 

A server scheduling policy at each time slot should decide how to allocate servers from set 
/C to the queues in set A/". This must be accomplished based on the available information about 
the channel state of the system at time slot t (i.e., C{t)) and also the queue length state at the 
beginning of time slot t (i.e., X{t — 1)). Therefore, at the beginning of each time slot t the 

^We assume that at each time slot departures occur first and then at the end of the time slot new arrivals are added to the 
queues. 
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scheduler has to determine an allocation matrix I{t) G X where 

I = Unxk = {Ium) , InM G {0, 1} I ^-^^ < 1 Vfc G /C 

I n=l 

and X is called allocation matrix space. We observe that each matrix in X can have at most a 
single "1" in each of its columns. The queue length vector evolves with time according to the 
following rule. 

= - 1) - {C{t) ® I{t)) il)^ + A'it) (2) 

where (•)+ is defined as follows: For an arbitrary vector v of size |^|, (v)^ is a vector of the 
same size whose i'th element, {v)^, is defined by 

f if < 

{v)t={ • (3) 

y Vi if > 

IV. Background 

A. Strong Stability Definition 

We now introduce the definition of strong stability for a queueing system [[T|, [[2). Other 



definitions can be found in pO] , ^9j, pT| , p2] . Consider a discrete time single queue system 
with an arrival process A{t) and service process fi{t). As we mentioned earlier, the arrivals are 
added to the system at the end of each time slot. We can see that the queue length process X{t) 
evolves with time according to the following recursion; 

X{t) = {X{t-l)-^{t)y + A{t) (4) 

Definition 1: A queue evolving with time according to (|4]) is said to be strongly stable \2j if 

1 '"^ 

Iimsup-^£;[X(r)] < oc. (5) 



T=0 



Definition 2: A queueing network is said to be strongly stable \^ if all the queues in the 
system are strongly stable. 

In this paper, we will employ the definition of strong stability and in what follows we use the 
terms "stability" and "strong stability" interchangeably. In [[T|, Q, it was proved that if a queue 
is strongly stable and if for all t either E[A{t)] < A or E[iJi{t) - A{t)\ < D where A and D are 
finite given non-negative constants, then 

lim -E[X{t)] = 0. (6) 

t^oo t 
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B. Lyapunov Drift 

A very important and useful mathematical tool used in network stability analysis and stochastic 
control/optimization of wireless networks is the Lyapunov Drift technique introduced in [[T|, Q, 
pol l, p2l |, [[23]]. The main idea behind the Lyapunov stability method is to define a nonnegative 
function of queue backlogs which can be seen as a measure of the total aggregated backlog 
in the system at time slot t. Then, we evaluate the "drift" of such function in two successive 
time slots by taking the effect of control decisions (scheduling or resource allocation policy) 
into account. If the expected value of the drift is negative as the backlog goes beyond a fixed 
threshold, then the system is stable. This method was used in [[T|, Q, [[20), [|22|, [ [23]] to prove 
the stability of a number of queueing systems. 

For a queueing system with queues and queue length vector X{t) = ^^^(t)), 
the following quadratic function has been usually used in literature ( [[T|, [[2|, [ |20| |, [[22), [ [23] ]) 
as a Lyapunov function. 

N 

L{X) = Y,Xl{t) (7) 

n=l 

Assume that £'[X^(0)] < oc, Vn G A/" and X{t) evolves with some probabilistic law (not 
necessarily Markovian). Then, the following holds [[2). 

Lemma 1: If there exist constants B > Q and 5 > such that for all time slots t we have 

N 

E[L{X{t + l))-L{X{t))\X{t)]<B-5Y,X^{t), (8) 

n=l 

then the system is strongly stable and further we have 

t-l N 

limsup-VVE[X„(r)]<-. (9) 

The proof of the lemma can be found in [|T|, [|2). The left hand side of expression ([8]) is usually 
called Lyapunov drift function which is a measure of the expected value of changes in the backlog 
in two successive time slots. From Lemma [TJ we can easily see the idea behind Lyapunov method 
in stabilizing queueing systems. It is not hard to show that, when the aggregated backlog in the 
system goes beyond the bound y , the Lyapunov drift in the left hand side of ([s]) will be negative, 
meaning that the system receives a negative drift on the expected aggregated backlog in two 
successive time slots. In other words, the system tends toward lower backlogs and this results 
in its stability. 
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C. Fundamental Concepts of Polytopes 

We present a brief review on convex polytopes and fundamental properties of them. These 
concepts will be needed for specifying the stability region of MQMS system in section [V} 



Definition 3: A convex polytope is defined as the convex hull of a finite set of points [24|, 



According to Weyl's Theorem [24|, a polytope in always can be expressed by a set 
7^ = {x G I a^x^ < /3i for £ = 1, 2, L} for some positive integer L and ai G and 
/5^gM. 

Definition 4: The dimension of a polytope V is represented by dim(7^) and is defined as one 



less than the maximum number of affinely independent points in V p4] , p5] |. 
A polytope V C is said to be full dimensional if dim{V) = A^. 

Dimension Theorem: For a polytope V C M^, dimension of 7^ is equal to minus the 
maximum number of linearly independent equations satisfied by all the points in V. 

Definition 5: For given a and b, equality ax^ < b is called valid for polytope V if for every 
point Xq G V, axQ < b. 

Definition 6: A face of polytope V is defined as J' = {x G P | ax^ = b} where inequality 
ax^ < 6 is a valid inequality for V. 

We call the valid inequality axj = b a face defining hyperplane for V if its associated face 
is not empty. Therefore, ax^ = 6 is a face defining hyperplane for V if it intersects with V at 
least at one point. Note that V has finitely many faces. However, the face defining hyperplanes 
of a polytope can be infinite. 

Definition 7: A facet of polytope 7^ is a maximal face distinct^ from V [25]. Hence, all faces 
of V with dimension dim (7^) — 1 are called facets of V. 

For a polytope 7^ = {x G | a^x^ < (5^ for £ = 1,2,...,L} an inequality is redundant if 
polytope V remains unchanged by removing the inequality. 

Redundancy Theorem in Polytopes [24]: Face defining hyperplanes describing faces of dimen- 
sion less than dim (7^) — 1 are redundant. 

Redundancy theorem states that to describe a polytope completely, only facet defining hyper- 
planes are sufficient. 

^Maximal relative to inclusion 
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V. Linear Algebraic Representation of Stability Region Ploytope 

In this section, we first introduce the convex hull representation of the stability region [[T|, Q. 
Then, we derive the necessary and sufficient conditions for the stability of MQMS system with 
stationary channel processes and general arrival processes with finite first and second moments. It 
will be shown that for stationary packet arrival processes, these conditions establish the stability 
region polytope. We determine a linear algebraic representation of the stability region in which 
we explicitly determine the coefficients of the linear inequalities describing the stability region. 

A. Stability Region Geometry 

Consider the class of deterministic policies Q where at each state of the system, each policy 
g ^ Q allocates the servers according to a predetermined allocation matrix depending on the 
channel state matrix (i.e., matrix C{t)). More specifically, for each policy g ^ there exists a 
one to one mapping from the channel state space S to the allocation matrix space X, namely 
/(^) : S I — y I. Therefore, each deterministic policy g is specified by \S\ allocation matrices 
s e S. The scheduler observes the state of the system and then based on the observed state 
5, it allocates the servers according to an allocation matrix Is^\ Note that set X is a finite set 
and since the channel state space is also finite, set Q is finite. 

Each deterministic policy g provides an average transmission rate for each queue n. Let Rn"^ 
denote the time averaged transmission rate provided to queue n and i?^^) = {R[^\ i?2^\ R^^^) 
the vector of average transmission rates. By conditioning on the channel state of the system, it 
is not hard to see that for each deterministic server allocation policy g, we have 

^^'^= {Y.^s{Cs®li'^)ll^ . (10) 

Note that each rate vector i?^^) determines a single point in M^. Now, consider the convex 
hull of all the points R^^)^ g e Q in M^, i.e., 

n = conv.hull (11) 

Each point in TZ, let say R^ can be represented by a convex combination of a finite set of 
points, R^^\g G i.e., 

\G\ 1^1 

i?" = J2p^R'^ .g^eg , = 1. > o. (12) 

i=i 1=1 
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Hence, IZ contains all the achievable transmission rate vectors of MQMS system. To achieve the 
transmission rate vector W ^IZ ii is enough to select policy gi with probability pi, i.e., in pi 
fraction of time slots. In other words, all the transmission rate vectors eTZ are achievable by 
applying a randomized policy that at each time slot selects policy gi with probability pi. Thus, 
the following Lemma follows. 

Lemma 2: The set of achievable transmission rate vectors TZ is specified by a polytope. 
Proof: The lemma follows directly from the definition of polytope in ([3]), equation ( [TT] ) and 
the fact that set Q is finite. ■ 

We will denote the achievable transmission rate polytope by V. According to the definition of 
stability region and the discussion in [[T|, Q, [ [TT| , (20) regarding the network stability region, 
we can conclude that polytope V specifies the stability region of MQMS queueing system with 
stationary channel distribution and stationary arrival processes (in this case E[An{t)\ = and 
A = (Ai, A2, Aat)). Specifically, it has been shown that if the system is stable, we should 
have A G P and also if A G P — hound{V), then there exists a server allocation policy that can 
stabilize the system. 

Although stability region of MQMS system is described by ( [T2| ), it cannot be applied as 
the constraints of utility optimization problem ([T]). In fact, ( [T2| ) only provides us an implicit 
description of stability region. To solve network optimization problems like ([T]), we should 
characterize the stability region by a set of linear/non-linear convex inequalities (or equalities) 
which is given in the following. 

B. Necessary Conditions for the Stability of MQMS System 

We first give an outline of the steps we take to find the necessary conditions for the stability 
of the system. 

• Using Lemmas |3]|5| we will characterize all the face defining hyperplanes of the stability 
region polytope. More specifically, we will determine the normal vector associated to each 
hyperplane as well as a single point where the hyperplane touches the stability region 
polytope. We show that each vector in is associated to a face defining hyperplane. 

• The stability region is a polytope and therefore can be characterized by the finite set of 
its facet defining hyperplanes. We will introduce a subset of M^, namely V C and in 
Lemma [6] we show that the vectors outside set V cannot be associated to a facet defining 
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hyperplanes of the stability region polytope. Note that set V is an infinite set itself. 

• In Lemma |7} we will show that although set V is an infinite set, it is finite up to multiplication 
of vectors by positive scalars, i.e., all the vectors in V can be produced by multiplying a 
positive scalar to a vector in a finite set V . 

• Using Lemmas [3] to |7} we will prove Theorem [T] which states the necessary conditions for 
the stability of the system. We also argue that the derivation of set V becomes difficult for 
large and M. We introduce a finite superset of set V whose elements can be specified 
easily. 

We introduce the departure matrix HNxK{t) = {Hn^k{t)),n G A/", G /C in which H^^kit) 
(Hn,k{t) < Cn^k{t)) represents the total number of packets served by server k from queue n at 
time slot t. Thus, the departure process from queue n at time slot t would be Hn^kit)- The 

following equality illustrates the evolution of queue length process with time. 



K 



X^{t) = X^{t - 1) - V H^Mt) + Mt) ^neAf (13) 



k=l 



For a strongly stable MQMS queueing system we can prove the following Lemma. 
Lemma 3: If the MQMS system is strongly stable, then for any vector a = (ai, 0^2, <^Ar) ^ 
we have 



aE[A{t)]^ = aE[H{t)]rj^. (14) 



Proof: The recursion ( [13] ) will result into 

t K t 

X^{t) = X.(0) -J2Y.^r^Ar) + (15) 

T=l k = l T=l 

By multiplying vector a to the queue length vector X{t) we have 

t t 

aX^{t) = aX^{0) - ^ aH{r)ll + ^ aA^{r). (16) 

T=l T=l 

Taking the expectation from both sides, dividing by t and then taking the limit as t goes to 
infinity, we will have the following. 

lim -aE[X{t)]^ = lim -aE[X{0)]^ - lim -S" aE[H{T)]ll + lim -S" aE[A{T)]^ (17) 

r=l r=l 

According to ^ and the assumption that E[X{0)] < oc, the left hand side term and the first 
term in the right hand side term are equal to zero and therefore the result follows. ■ 
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A direct result of Lemma [3] for a single queue n is the following. 



E[Ar,{t)] = E 



K 



.k=l 



(18) 



In other words, for a strongly stable queue n, the time averaged total expected arrival to queue 
n is equal to the time averaged total expected departure from queue n. 



Lemma 4: If the MQMS system is strongly stable, then E[A{t)] G V. 
Proof: The proof follows directly by considering ( [T8] ) and the fact that 



K 



k=l 



<Rn Vn G A/". 



(19) 



Based on Lemmas [3] and |4} we can prove the following Lemma. 
Lemma 5: If the MQMS system is strongly stable, then for all a G 

aE[A{t)Y <^^s max ® /) 1^) 



(20) 



Proof: Consider a G M^. Since the system is strongly stable, from Lemma [s] we have 



lim - ^EiMr)r = lim y ^ a^[//(r)] 1 



T=l T=l 

By conditioning on the channel state process, we will have 

t 



1^ 



limi^«^[i/(r)]l 

r=l 

= /i^ i E E ^^"^[^(^) I ^(^) = 

r=l sG<S 

^ /i^ i E E ^^(")^^ ® ^(^)] 

r=l s65 
1 * 

< lim - E E 7¥ ((")^(^^ ® 

r=l sG<S 

= J]7r,inax((a)+(C,®/)1^) 



sG5 



Now, consider two different cases: 



a G M^^: The result follows directly from (21) and (22) since = a. 



(21) 



(22) 
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a ^ M^: In this case, from (21) and ( [22| ) we have the following inequality. 

aE[A{t)Y <^7rs max ® 

However, since (a)+ G and according to the previous case we also have 



(23) 



(24) 



Noting (24) and the fact that E[A{t)] G M^, we conclude that (23) is a redundant inequality. 



Since a G M^, the set of inequalities in (20) forms an infinite set. Each inequality in (20) 



determines a valid inequality for poly tope V. 



Fact: The hyperplanes associated to the valid inequalities of ( [20| ) are all face defining hyper- 
planes of poly tope V. 

To show this fact, let denote a set of allocation matrices {I^.s G S} that maximize the 



right hand side of (20), i.e.. 



(25) 



is not unique and there may be more than one set of allocation matrices of whose elements 
maximize a{Cs®I)lJc. According to (10) and 

111])' iT.ses {Cs ® I^) llf G V. On the other 
hand, since Y.^^^ 7Tsa{Cs®Is)lK = ^ T^ses ^s{Cs®I^)Ik, point Y^ses ^s{Cs®Is)Ik is located 
on the hyperplane associated to ( [20| ). Therefore, the set of inequalities in ( [20| ) determines all 
the non-empty faces of polytope V. To clarify what we proved in Lemmas [3] to |5] consider the 
following example. 

Example: Consider a system with = 2, M = 1 (ON-OFF channels) and K = 1 (see 
Figure |2]). Assume that queues 1 and 2 are connected to the server with probabilities pi and p2, 



respectively. This system with N >2 queues was studied in p9] when studying dynamic server 
allocation to a set of parallel queues. Here we consider a system with just 2 queues so that we 
can illustrate the stability region in two dimensions. For such a system it was proven that the 
stability region is described by the following set of inequalities. 

Ai < pi 

A2 < P2 

Ai + A2 < pi + P2 - P1P2 (26) 
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Ai 



As 




Fig. 2: Two-queue single- server example 



In Lemmas [3] to |5} we characterized the stability region by an infinite number of its face 
defining hyperplanes. Figure |3] depicts the stability region given in ([26]) as well as some of its 



face defining hyperplanes derived according to (20) for some sample values of a = (ai, 0^2) 




Fig. 3: Stability region for the queueing system of Figure j2j 



According to the redundancy theorem in poly topes (refer to Section |IV-C[ ), not all the face 
defining hyperplanes of a polytope are required to describe a polytope. In particular, just the 
inequalities corresponding to the facets of a polytope are sufficient to characterize a polytope. In 
the following, we will characterize the facets of polytope V. The analysis follows and constitutes 
the main contribution of the paper. 

Let V denote the set of vectors a G with the following property, i.e.. 



(27) 



A vector a belongs to set V if for any partitioning of elements of vector a into two non-empty 
disjoint (sub)vectors in which all the elements of each vector are not identically equal to zero, 
there exists at least one non-zero element in each (sub)vector, the ratio of which is equal to 
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the ratio of two non-zero elements of set M. To clarify the definition of set V consider the 
following example. 

Example: Let M = {0, 1, 2} and = 4. Now, consider vector ai = (1, 2, 5, 10). We can 
partition ai to ai^^ = (1^ 2) and o^i^g^^ = (5, 10). Note that there exist no two elements one in 
ai^^ 2} the other in ^1^3 whose ratio is 1 or 2 and therefore, ai ^ V. Another example is 
vector a2 = (0, 2.5, 0, 0) in which for any partitioning of the vector into two non-empty disjoint 
(sub)vectors one will always have all the elements identically equal to zero. Therefore, 0^2 G V. 

In the following, we will show that any vector in — V cannot form a facet defining 
hyperplane for polytope V and therefore the inequalities derived by the vectors of set — V 
are redundant. This is shown in the following lemma. 

Lemma 6: If a G — V, the hyperplane associated to the valid inequality of (20) is not a 
facet defining hyperplane of V . 
The proof is brought in Appendix |A} 

Note that V is an infinite set. However, we can show that set V is finite up to multiplication of 
vectors by positive scalars. To show this, we define set >V = |z G Z+ | z = Y[J=i ^ -^j 

and = {(^1, 0^2, g^n) \ an G W}. Then, we can prove the following Lemma. 

Lemma 7: There exists V C such that any vector a e V can be written as a = g/3 for 
some vector /3 and scalar q > 0. 
The proof is given in Appendix [B} 

Recalling the two-queue single- server example of Figure [2} using Lemmas |6] and |7] we 
can characterize a finite subset of that can produce all the normal vectors of the facet 
defining hyperplanes of the stability region polytope. For the aforementioned example since 
M = 1 we have W = {0,1} and therefore = {(0, 0), (0, 1), (1, 0), (1, 1)} and therefore 
V C {(0,1),(1,0),(1,1)}. For this simple example since M = 1 (as we will see in Corollary 
[2] below), we have V = — {0^}. The facet defining hyperplanes of the stability region for 
this example as described by the vectors in the set V are shown with bold red lines in Figure [s] 

Fact: The size of >V is |W| = {^^^~^) + 1. 

In order to prove this fact, note that M has M non-zero elements. Each non-zero element 
of W comes from the multiplication of — 1 non-zero elements of M. This is equivalent 
to the counting problem of choosing — 1 balls from M distinctly marked balls without 
ordering and with replacement which is equal to |>V| = [26]. Since G W we 
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must increase this number by 1. Hence, = {{^^^^^) + l)^. Since V C - {0^}, 

therefore < |>V^| — 1 (excluding 0^ as it results into the obvious equality = 0). 

In the theorem that follows, we introduce the necessary conditions for the stability of MQMS 
system with stationary channel distribution. 

Theorem 1: If there exists a server allocation policy under which the system is stable, then 



(28) 



Proof: The proof follows directly from Lemmas [5] to |7| In fact, from Lemma [5] we can 
characterize polytope V by an infinite number of inequalities (a G M^). In Lemmas |6] and [t] 
we showed that not all of such a vectors are necessary to describe the facets of polytope V. 
More specifically, we proved that just a G 1/ (|1/| < oc) are sufficient to describe the facets of 
polytope V and other a vectors make redundant inequalities and therefore the theorem follows. 



According to Theorem [TJ in order to characterize polytope V, we have to specify all the 
elements of set V which can be obtained after considering all the possible vectors a G W^— {0^} 



and removing redundancies following ( [27] ). This can be achieved numerically. We may also avoid 
the numerical computation to derive set V and instead use set — {0^} which is a finite 
superset of V that contains redundancies (produces redundant inequalities). Although by using 
— {Ojsi} we may obtain some redundant inequalities, since |>V^ — {0^}| < oc we still have 
a finite number of inequalities to describe polytope V. Table depicts t/ and |>V^| - 1 for some 
sample cases TV = 2, 3 and M = 1, 2, 3, 4. 

TABLE I: \V\ and |>V^| - 1 for = 2, 3 and M = 1, 2, 3, 4 





M = 1 


M = 2 


M = 3 


M = 4 


|V^|,|W^|-1 


iV = 2 


3,3 


5,8 


9,15 


12,24 


N = 3 


7,7 


25,63 


109,342 


253, 1330 



Corollary 1: For an MQMS system with stationary arrival processes, we have E[A{t)] = A 
(Ai, A2, Aat) and the necessary conditions for the stability of the system would be 



(29) 
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Corollary 2: For an MQMS system with ON-OFF channels, we have W = {0, 1} and 
therefore the necessary conditions for the stabihty of the system would be 

aE[A{t)Y <^^s max ® , a G {0, 1}^ - {0^}. (30) 

In a special case considered in where the channels are modeled by independent Bernoulli 
random variables with E[Cn,k{t)] = Pn,k^ the necessary conditions for the stability of the system 
are given by 

K 

Y^E[A^{t)]<K-Y,\{{'^-Pn.k) 'iQ<^M. (31) 

n^Q k=l n^Q 

In this case, the total number of inequalities needed to describe polytope V is equal to 2^ — 1. 



C. Sufficient Conditions for the Stability of MQMS System 

Consider a server allocation policy which determines the allocation matrix at each time slot 
t by solving the following maximization problem. 

I{t) = argmaxX(t - ® (32) 

This policy is called Maximum Weight (MW) and was introduced in [[T0|, [[T3|. According to the 
constraints on allocation matrix /, we can easily conclude that MW policy allocates the servers 
by the following rule: At each time slot t, each server k is allocated to the queue n that achieves 
the maximum Xn{t — l)Cn^k{t). 

In the special case where the channels are ON-OFF, the policy is the same as AS/LCQ (Any 



Server/Longest Connected Queue) introduced in p8] |. In the following, we derive the sufficient 
conditions for the stability of our model and prove that MW stabilizes the system as long as 
condition ([33]) below is satisfied. We also derive an upper bound for the time averaged expected 
number of packets in the system. 

Theorem 2: The MQMS system is stable under MW policy if for all t 

aE[A{t)]^<J2^srnd.x{a{Cs®I)ll), aeV. (33) 
Furthermore, the following bound for the average expected "aggregate" occupancy holds. 

,^..,,l^±ElX„ir)l<!l3.^pM^ (34) 

r=0 n=l 
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In (34), 5 is the maximum positive number such that /or all t we have E[A{t)\ + 51^ G P. 
The detailed proof is given in Appendix [C} 

Corollary 3: For an MQMS system with stationary arrival processes, we have E[A{t)] = 
A and the sufficient conditions for the stability of the system under MW policy are aX^ < 

According to Corollaries [T] and |3] and the definition of network stability region, we can conclude 
that for an MQMS system with stationary channel distribution and stationary arrival processes. 



the stability region is characterized by ( [29| ). The stability region described in ([29]) is a poly tope 
which can be imagined for 2 and 3 dimensional systems, i.e., = 2,3. Figure |4] shows the 
stability region for = 2, 3 and M = 1, 2, 3, 4. In all the cases K = 3, 

Consider an MQMS system with ON-OFF channels. For such a system, MW policy is 
equivalent to AS/LCQ policy introduced in JTS] ]. AS/LCQ policy takes an arbitrary ordering 
of servers (to be allocated to the queues) and then for each server allocates it to its longest 
connected queue (LCQ). Note that in such a system, W = {0, 1} and therefore we have the 
following corollary. 

Corollary 4: The MQMS system with ON-OFF channels is stable under AS/LCQ if for all t 

aE[A{t)]^ < max {a{Cs ® , a e {0, 1}^ - {0^}. (35) 

In a special case where the channels are modeled by independent Bernoulli random variables 
with E[Cn^k{t)] = Pn^k^ AS/LCQ stabilizes the system if for all t 

K 

Y,E[A^{t)\<K-Y,\[{l-p^,k) yQ^Af. (36) 

neQ k=l neQ 

D. Numerical Example 

In this section, we will consider a simple numerical example of MQMS system with N = 2, 
K = 3 and M = 3. For such as system the size of channel state space is 4^ = 4096. We have 
chosen the channel distribution randomly. We also assume that the utility functions associated 
to queues 1 and 2 are log(l + lOri) and 10r2, respectively^. Therefore the total utility function 
is /(ri, = log(l + lOri) + 10r2. We also assume that the packet arrival processes to queues 
1 and 2 are Poisson distributed with rate 5 packets/time slot. First, we are going to characterize 



^Linear and logarithmic utility functions were studied in |1|, ||2j 
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(a) iV = 2, M = 1 




(c) = 2, M = 2 



2 2.5 3 3.5 



(e) AT = 2, M = 3 
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(b) iV = 3, M = 1 



(d) iV = 3, M = 2 
















C [. 1, ■ 







(f) AT = 3, M = 3 




(g) AT = 2, M = 4 



(h) AT = 3, M = 4 



Fig. 4: Stability region for TV = 2, 3, M = 1, 2, 3, 4 and X = 3 
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the stability region for this example using the theorems given in the previous subsections. Then, 
we will correlate our findings with the results from [[T|, Q, Q. 

For the described system we can easily check that W = {0, 1, 2, 3}. Thus, >V^ — {(0, 0)} = 
{(0, 1), (0, 2), (0, 3), (1, 0), (1, 1), (1, 2), (1, 3), (2, 0), (2, 1), (2, 2), (2, 3), (3, 0), (3, 1), (3, 2), (3, 3)} 
whose size is 15 (as mentioned in Table |I]). Note that some of these vectors make redundant in- 
equalities. For example using (1, 1), (2, 2) and (3, 3) as the vector a in ( [29] ) will result in the same 
inequalities. By removing the redundant vectors we obtain set V = {(0, 1), (1, 0), (1, 1), (1, 2), 



(1, 3), (2, 1), (2, 3), (3, 1), (3, 2)}. Using (29), the stability region is characterized by the following 



set of inequalities which is also shown in Figure [5j 

ri < 4.4792 n + 2r2 < 10.2893 2ri + 3r2 < 16.4564 

r2 < 4.4912 ri + 3r2 < 14.6002 3ri + r2 < 14.5803 

ri + r2 < 6.3577 2ri + r2 < 10.2874 3ri + 2r2 < 16.4639 (37) 
Therefore, the utility optimization problem is specified as 

Maximize : log(l + lOri) + 10r2 (38) 



Subject to : All the constraints in (37) 



< ri < 5 , < r2 < 5. 

The solution to the above problem is = (r^, = (1.1266, 4.4912) and is indicated by the star 
in bold in Figure [5j As we expect the optimal point is located on the boundary of the stability 
region. For each queue, we can use a leaky bucket and adjust the admitted rate to each queue 
to the optimal value we obtained from the optimization problem. 

We have also simulated the flow control strategy CLC2b (Cross-Layer Control 2-b) which was 
introduced in [[T|, Q, [[6| (Due to lack of space the reader is suggested to refer to the literature 
regarding this algorithm). CLC2b was proposed as a flow control strategy for a general queueing 
system for which the stability region is unknown. The algorithm assumes that An{t) < K^^^ for 
all t. Assume that the backlog in the transport layer of queue n at the beginning of time slot 
t is denoted by Ln{t). CLC2b algorithm is a dynamic flow control algorithm in the transport 
layer in which at each time slot it determines how many packets to admit to queue n (in the 
network layer) from the transport layer (i.e., from L^(t) packets in the transport layer waiting 
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1 - 
0.5 - 



0.5 1 1.5 2 2.5 3 3.5 4 4.5 

n (packets/time slot) 

Fig. 5: Stability region for = 2, X = 3, M = 3 with a random channel distribution (the 
optimal solution and the solution of CLC2b algorithm are also depicted in this figure) 

for admission). Assume that this number is denoted by r^(t). Therefore, the goal of CLC2b 
algorithm is to determine r^(t) such that /(f) is maximized while the queues are kept stable. 
CLC2b algorithm incorporates two novel notions, namely auxiliary variables and virtual cost 
queues. In the following, we will briefly review them. 

The utility optimization in ([T]) can be transformed into the following stochastic optimization 
problem by introducing the auxiliary variables 7^, n G A/" (ij. @? [[6|. 

N 

Maximize: ^^fniln) (39) 

^'^ n=l 

Subject to: ^ 7n 

r = (ri,r2, -...Vn) G A 
0<r^< Vn = 1,2, ...,A. 

For each queue n a virtual cost queue with virtual arrival process 7n(t) and virtual service 
process r^(t) is defined (See Figure [6] and refer to [2]). Note that the service process of the 
virtual queue n is exactly equal to the admitted packet process to the actual queue n. Assuming 
that the length of this virtual queue is denoted by Yn{t) at time slot t, the following equation 
reveals the evolution of this process by time. 

Yn{t) = {Y^{t - 1) - rn{t)y + 7n(t) (40) 
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An{t) — 



Yn{t) 



Transport Layer 




Fig. 6: CLC2b algorithm framework 

It was shown in [fT|, Q that if a control strategy stabilizes both the actual and virtual queues 
then the resulting averages and 7^ will satisfy the constraints of the optimization problem 
([39]). The flow control strategy of CLC2b is specified as follows. 

• Every timeslot and for each queue n observe Xn{t — 1) and Yn{t — 1) and determine 

min{L,(t) + if vYnit - 1) > - 1) ^^^^ 

otherwise 
Choose jn{t) (the input process of queue n's virtual queue) as the solution of the following 



Tnit) = 



Optimization problem and update Yn{t) following (40). 



Maximize : 

7 



(42) 



VfM-VYn{t-l)^ 
Subject to: < 7 < i?;^'^^. 

• Use the MW policy for resource (server) allocation and update 
fn{-) is the utility function associated to queue n and is a positive weight that satisfies < 
77 < 1. It determines the relative weight of the virtual queue in stabilizing the system. V is 
a control parameter that affects the proximity of (the solution of CLC2b algorithm) to the 
optimal point r~^. It was shown that CLC2b can stabilize both the actual and virtual queues while 



guaranteeing a lower bound for the achieved utility function (given by ([43])) and also an upper 
bound for the system backlog in the network layer (given by (|44|)). 

D 



liminf /(r) > /(r^) - 



^ t N 

limsup-^^X^(r) < 



D + VG 



(43) 



(44) 



r=0 n=l 

Where D, G and /x are positive constants and are dependent to the statistical properties of the 
system. We can observe that as the control parameter V increases the difference between /(f) 
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and the optimal utility /(r^) decreases as ^. However, the closer /(f) and f{r~^) are the larger 
the backlog is in the system! 

We have simulated the CLC2b algorithm with r] = 1. Furthermore, we performed a simulation 
where the admitted rates to the queues are controlled by leaky buckets outputting the optimal 
rates = (1.1266, 4.4912) as determined by our analysis earlier. We compared the two systems 
by computing the following measures versus V. 



• Difference of the average queue occupancy of the systems, i.e., XcLC2b — Xopt shown in 



Figure 7a 



• Percentage of the normalized difference of the utility function for r and T , I.e., ^ ^ ^ ^ 

100 shown in Figure |7b} 
We can see from the graphs that the queue occupancy under CLC2b algorithm is growing with V 
linearly (as expected by Q) while the utility associated to the CLC2b algorithm is converging 
to the optimal point like ^. The solution of CLC2b algorithm for different V parameters are 
also shown in Figure [5] (shown by the various asterisks). We can observe that as V gets larger 
the solution of CLC2b algorithm converges to the optimal point. 




(a) Difference of the average queue lengths for the optimal (b) Percentage of the normalized difference of the utility 
solution and the CLC2b algorithm function for r and r'^ 



Fig. 7: Comparison of the performances of CLC2b algorithm with the optimal solution 

E. Stability Region for Fluid Model MQMS Systems with Stationary Continuous Channel Dis- 
tribution 

We will consider a time slotted fluid model MQMS system with stationary channel distribution. 
In this case, the amount of work that arrives into (and departs from) the queues is considered 
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to be a continuous process. We also assume that the channel state of the link from each queue 
to each server is modeled by a continuous random variable. The channel state matrix in the 
fluid model MQMS system is defined as C{t) = {Cn,k{t)) ,n e JV,k e ]C,Cn,k{t) G M+. We 
assume that the channel process follows a stationary distribution fc{t){c) where c G M^^^; 
i.e., /c(-) is the joint distribution of all Cn,k{t) variables. We can easily check that in the fluid 
model MQMS system, set V is equivalent to as in the fluid model MQMS system, set M 
is replaced with M+. Therefore, stability region for the fluid model system is characterized by 
the following set of inequalities. 

POO POO POO 

aX^ < / • • • / max {a{c ® I)!^) 

xfc{t){c) dci^i • • • dcn,k-i dcn^k G (45) 
Note that in order to characterize the stability region of the fluid model system we need 



to compute an infinite number of nested integrals in ( [451 ). f^^^' the stability region of fluid 
model system is characterized by an infinite number of half spaces, hence the stability region 
is a convex surface. In this case, depending on the channel distribution and dimension of the 
system, we may characterize the stability region by a finite number of non-linear inequalities 
instead of infinite number of linear inequalities as we show in the following example. 

Example: Consider a fluid model MQMS system with two queues and one server. Assume 
that the channel state variables Ci^i(t) and C2,i(t) are independent and follow exponential distri- 
bution with means jii and /X2, respectively. Such a model may be used for slow Rayleigh fading 
channels in low SNR regimes where SNR follows exponential distribution and the approximation 



log(l+x) X is used for small positive x. According to (45), the stability region is characterized 
by 

^ POO POO C2,i 

aiAi + ^2^2 < / / maxjaiCi^i, a2C2,i}e e ^2 dc2^idci^i Vai,a2 gM+ (46) 

M1M2 Jci,i=0 Jc2,i=0 

We can write the right hand side of (|46]) as 



1 r^'^i'i _fM ^2,1 

— / / <^iCi,ie ^1 e ^2 dc2^idci^i 

lM2 Jci,i=0 Jc2,i=0 

2 rOO POO C2,i 

H / / <^2C2,ie ^1 e ^2 rfc2,irfci^i 

M1M2 Jci,i=0 Jc2,i = ^ci,i 

= / Ci^ie ^1 / e ^2 dc2^idci^i 

M1M2 Jci 1=0 ' Jc2 1=0 
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A2 < (yr^) (2 - 



/ii = 2 and /i2 = 1 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Fig. 8: Stability region for /ii = 2 and /i2 = 1 



+ 



Q^2 



e Ml 



Jci,i=0 

( Ml - 



C2,l = ^Cl,l 



+ Q^2 



C2,ie ^2 rfC2,lrfCi,i 

/ aia2/X2/xi 



+ 



0^2/^2 



Therefore, all the ordered pairs (Ai, A2) = f/xi - , "'j^^ , ^i^2Mim + "^^2 



(47) 



char- 



acterize the boundary of the stability region. However, we can write A2 based on Ai, fii and /i2 
as follows. 



Ml 



(o^i/ii + ^2^2)^ 



So, 



Q^2M2 



'1-^ and 



Ai 




Therefore, 



A2 = M2 Wl 




(48) 



As we can see from equation ( [48] ), the stability region in this example is characterized by 
just one non-linear inequality A2 < /X2 (^/^^^^^) ~ -^1 — and two linear inequalities 
Ai > and A2 > 0. The stability region for this example for /xi = 2 and //2 = 1 is illustrated in 
Figure [8| The characterization of the general fluid model stability region is beyond the scope of 
the paper and left as a possible future research problem. 
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VI. Conclusions and Future Work 

In this paper, we introduced a linear algebraic representation of the network stability region 
(capacity region) polytope of multi-queue multi-server (MQMS) queueing system with stationary 
channel distribution and stationary arrival processes. To this end, we obtained the necessary and 
sufficient conditions for the stability of the system for a general arrival process with finite 
first and second moments. For stationary arrival processes, we showed that these conditions 
establish the network stability region of the system given by ([29]). For the stability region 
polytope, we explicitly determined all the coefficients a e V of all the half spaces which 
are required to characterize the stability region polytope. We also argued that in general it may 
be (computationally) hard to quantify set V. In this case, although we may add some redundant 
inequalities, we can use a superset of V, namely the set — {0^} (V C — {0^}) instead 



of V in ( |29| ). An upper bound was also obtained for the average queueing delay of Maximum 
Weight (MW) server allocation policy which is a throughput optimal policy for MQMS system. 
We finally considered the stability region for a fluid model MQMS. In this case, we determine 
the stability region by an infinite set of linear inequalities given by ( [45] ). By use of an example we 
showed that depending on the channel distribution and the number of queues, we may characterize 
the stability region by a finite set of non-linear inequalities instead of infinite number of linear 
inequalities. However, the general problem of stability region characterization for the fluid model 
MQMS can be considered as a possible future research problem. 

Appendix A 
Proof of Lemma [6] 

Proof: As we explained before, denotes a set of allocation matrices of {/^, 5 G 5} that 



maximize the right hand side of (20) and we discussed that is not a unique solution of 



and there may be more than one set of allocation matrices of whose elements maximize ( [25] ). 



Let = {Xf , 1 < i < 11^1} denote the set of all distinguished solutions of (25). Obviously 
< oc, since set I is finite. Note that each solution = G 5}, 1 < i < is 

corresponding to a deterministic policy with the average transmission rate vector 

R?=(^7i4C,®i:^,)ll^ (49) 
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and therefore, according to ( pTj ) each solution of ( |25| ) is associated with a vertex of polytope 
V. Since there are distinguished solution of ([25]), vertices of polytope V are on the 
hyperplane associated to half space ( [20| ), i.e., 

aE[A{t)Y = ^^s max {a{Cs ® /) 1^) (50) 



Therefore, each face of V defined by hyperplane (50) is represented by F^, = conv.hulli<i<\ic.\Rf. 
Since polytope V is full dimensional, according to Definition [7} each facet of V must be of 
dimension — 1. In the following, we will show that the dimension of face is less than 
N — 1 if a e — V and therefore the hyperplane (50) associated to such a is not a facet 
defining hyperplane. 

Consider a vector a G — V. For a we have the following property: 

y{iGU^jGW^ m, n G A^, Q^i, Qfj, m, n 7^ 0) aiiri aju (51) 

In other words, there exists a partitioning of vector a into two (sub)vectors au and auc such that 
no non-zero element of au is proportional to no non-zero element of auc by the ratio of two 
non-zero elements in M. Assume that Uq denotes the subset U satisfying ( [5T] ). Now, consider 



inequality ( |20| ). The right hand side of ([20]) can also be expressed as 

K 

TTs max {a{Cs ® = XI ((C'f )^ ® ^)> , (52) 

seS seS k=i 

where is the fc'th column of matrix Cg and is the set of all binary vectors of size with 

J:={J= ( Ji, J2, ^iv) I Jn e {0, 1} Vn G A/", (J, 1^) = 1}. (53) 

Each Js{a) = argmaxj^j^ (a, ((Cj-^)^ ® J)) is corresponding to the fc'th column of an in 
some G I^. Note that maxj^j (a, ((Cf^)^ ® J)) > 0. There are K\S\ of such maximization 
terms in the right hand side of ( [52] ). Each of the K\S\ maximization terms may have multiple 
solutions (non-unique solutions). Suppose that (o^) denotes the set of all distinct solutions to 

J^{a) = argmaxj^^ {a, {{Cj^f ® J)), i.e., 

i^si^) = Usi^) = argmax(a, ((C^f )^ ® ^)> G 5, G /C} (54) 
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For a given 5 G 5, k e JC with (a, ((Cf^)^ ® J^{a))) > 0, we can easily observe that for 
all the elements of ips{a) either we have XIwgz^o "^s^,ii(^) = 1 = 1- This result 

comes directly from pT] ) and states that for a given 5 G 5, G /C, /c^r ^2// the solutions 
of ^^^(a) = argmaxj^j (a, ((Cf^)^ ® J)) with maxj^j (a, ((Cf^)^ ® J)) > 0, the index 
associated to element "1" of Js{a) is either in Uq or Wq- other words, it is not possible 
to have the index of element "1" in Uq for some solutions and in Uq for other solutions of 
Jg{a) = argmaxj^j (a, {{Cj^)^ ® J)). Using this property, we can partition set of all ipg{a) 
with (a, ((Cj-^)^ ® Jg{a))) > into two disjoint subsets. We denote the set of all solution sets 
^l;^^{a) with (a, ((Cf )^ ® J^{a))) > and J^ueUo J^ui^) = 1 as X and set of all solution 
sets ^^^{a) with (a, ((Cf )T ® > and Y.u&a^^ J^ui^) = 1 ^s We introduce 

dimensional vectors a^^^ and a^c as follows. 





if neU 









a new 








Now, let us introduce the following two hyperplanes. 

a'^,E[Au{t)V = E vr,max(«,((Cfr® J)) (55) 

a'^.E[AAW= E vr,max(a,((Cff ® J)) (56) 



In the following, we will prove that all the vertices Rf, 1 < i < in (49) satisfy both (55) 



and ([56]). In other words, all the vertices Rf are located on both hyperplanes ( [55] ) and ([56]) and 



therefore on the intersection of both. The hyperplanes of ( [55] ) and ( [56] ) are in an dimensional 



space whose dimensions are at most — 1. Thus, dimension of their intersection is at most 
— 2. In other words, dimF^, = dim{conv.hulli<i<\io.\Rf) < N — 1. Therefore is not a 
facet for V. 

Note that columns of any allocation matrix are vectors Jg{a) G for all G /C. 

Consider a particular transmission rate vector Rf (equation ([49])) located on hyperplane ( [50] ). 



As we discussed before, each allocation matrix I^- in ( [49] ) is created by concatenation of a 
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set of vectors Jg^i{ct) G 'ipsi'^) the columns of Now, we can check that all the vertices 



Rf, 1 < i < |I"| in (|49|) satisfy both (|55|) and (|56|) as follows. 



and 



K 



seS k=i 



J2 TT, max (a, ((Cf)T® J)) (57) 



K 



= E^^E(<-((^'')^®^^(«)))= E ^.max<a,((C7fr®J)> (58) 

s^S k=l k,s-4^{a)^Bc 



Thus, all the vertex ( [49] ) satisfy both ( [55] ) and ( [56] ) and therefore their intersection. This shows 
that Fo, has dimension less than — 1 and therefore is not a facet of V. ■ 



Appendix B 
Proof of Lemma [7] 



Proof: Consider vector a = (ai, 0^2, <^Ar) ^ V. From (27), it is obvious that any vector 
derived by the multiplication of a positive scalar to a also belongs to V, i.e.. 



(59) 



Also note that vectors a and qa will result in the same inequalities in ( [20] ) and therefore the 
same face defining hyperplanes of poly tope V. Thus, we say that a and qa are equivalent and 
we write a = qa. 

Consider a vector a . Vector a may be a zero vector in which case it does not contributes 



to a face of V as the inequality ( [20] ) will result in an obvious equality = 0. Hence, we assume 
that a is not a zero vector. 

Now, we perform the following process on vector a. Pick a non-zero elements of a, let say a^^ 
and form a{e^} and o^A/'-{ei}- The size of these vectors are 1 and — 1, respectively. According 



to (27), the vector a^v'-jei} should be either a zero vector in which case we stop or there should 



be a non-zero element a^^ in a;v'-{ei} and non-zero mi,ni G M such that a, 



62 



ei 



mi 



. If SO, 



we will continue and form Q^{ei,e2} <^A/'-{ei,e2}- According to (27), the vector Q^A/'-{ei,e2} is 
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either a zero vector in which case we will stop or there should exist a non-zero a^^ G ajv'-{ei,e2} 



and non-zero m2,n2 such that either 



^6 CI ^2 ^3 C2 ^2 ^1 ni 722 

procedure until either for some < the vector Q^A/'-{ei,e2, -,ei.} ^ zero vector or we form 
vector Q^{ei,e2,. -,eiv} which case u = N.\n both cases, we can see that all the non-zero elements 
of a may be expressed by the multiplication of ai and the rational numbers derived by non-zero 
elements of set M, By dividing the vector a by scalar a^-^ and multiplying it by Hjli ^j' we 
will obtain an equivalent vector for a named (5 whose elements belong to the set of all integers 
which were derived by multiplication of — 1 elements of set M, i.e., 

a = (3 , G W = |z G Z+ I z = ]J m^-, m^- G Aij Vn G {1,2, ...,A^} (60) 

Therefore, each vector a G ^ is equivalent to a vector /5 whose elements come from set W. 
We define the set of all (5 vectors with the above property as V. Note that V C and is a 
finite set since W is finite. Thus, for each vector a there exists an equivalent vector (5 
and the result follows. ■ 



ap 



1712 



ap 



mi 1712 



We repeat this 



Appendix C 
Proof of Theorem [2] 

Proof: We will start with the Lyapunov function evaluation, we will use the quadratic 
function (|7]) as our Lyapunov function. The Lyapunov drift for two successive time slots has the 
following form. 



D{t + 1) = E[L{X{t + 1)) - L{X{t)) I X{t)] 



= E 



N 



Y,K{t + i)-xl{t)\x{i) 



n=l 



= E 



N 



Y,{Xn{t+l)-X^{t)f\X{t) 



n=l 



2E 



N 



J2XnmXn{t+l)-Xn{t))\X{t) 



n=l 



(61) 



For the the first term we have: 



E 



N 



Y,{Xn{t + l)-X^{t)f\X{t) 

N K 
= E Y,{Ar,{t+l)-Y,Hr.,,{t + l)f\X{t) 



n=l 



k=l 



= E 



N 



j2Ai{t+i)\x{t) 



n=l 
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' N K 






N 




2E 


n=l k=l 


+ 


E 


E 

n=l 


^^H^,kit + l)^ \Xit) 



N 



K 



n=l 

For the second term in ([61]) we have 

E 



(62) 



N 



n=l 



J2Xnit){X^{t+l)-X^{t))\X{t) 

' N / K 



E 



N 



n=l 



k=l 



N 



n=l 



n=l 



K 



Xn{t)Y,HnAt+^)\X{t) 



k=l 



(63) 



Therefore, the Lyapunov drift D{t + 1) can be bounded by 



N 



D{t^l)<NAl,, + Y,E 



n=l 



K 



N 



N 



+ 2 J] Xr,{t)E[An{t + 1)] - 2 J] ^ 



n=l 



n=l 



K 



Xn{t)Y,Hn^k{t+l)\X{t) 



k=l 



(64) 



In the following, we show that the Lyapunov drift in ( |64| ) is bounded as follows. 

2 



N 



D(t + i)<iVAL. + E^ 



n=l 



K 



Y.Cn^k{t+l)In^k{t+l)\ \X{t) 



Kk=l 



N 



N 



+ 2 Xn{t)E[An{t + 1)] - 2 J] E 



n=l 



n=l 



K 



Xn{t)Y,CnAt+^)InAt+^)\X{t) 



k=l 



(65) 



To prove ([65]), note that for each queue n e JV one of the following conditions is satisfied. 

• ^k=i Hn,k{t + 1) = J2k=i Cn,k{t + i)In,k{t + l)- this casc, wc cau easily check that for 
all n G A/" 



E[Xl{t + l)-Xl{t)\X{t)]<NA, 



2 

max 



+E 



K 



^k=i 



-2Xn{t)E[An{t + 1)] - 2E 



K 



X^{t)Y,Cn,k{t + l)InA'^+l)\X{t) 



k=l 



(66) 
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Eti HnAt + 1) < Ek=i CnAt + + 1): In this case, Eti + 1) = Mt) 

and + 1) = 74^(t + 1) as there are not enough packets in queue n to be server at time 
slot t + 1. It is not hard to check that inequahty (|66|) is also satisfied in this case. 



According to the above discussion, we can observe that inequality ([66]) is satisfied for all the 
queues and therefore ([65]) follows. 

K 

Using the the fact that ^ Cn,k{t + ^)In,k{t + 1) > 0, we get the following inequality. 



k=i 



N 



K 



N K 



5] I]C'n,fc(t + l)4,fc(t + l) < + + <{MKf (67) 



n=l \k=l 



^n=l k=l 



Hence, the Lyapunov drift ( [65] ) is bounded by 

D{t + 1) < NAl,, + {MKf + 2 Xn{t)E[Mt + 1)] 



N 



N 



n=l 



n=l 
K 



k=l 

By conditioning the last term of ( [68] ) on the channel state at time slot t + 1, we will have 

D{t + 1) < iVAL. + (MKf 



(68) 



-2E 



Xit)lA^{t + l)-J2^,iCs®I{t + l))ll] \Xit) 



(69) 



Note that the allocation matrix / in d69]) depends on the selected policy. According to (32), we 



can see that by selecting MW policy, the second term of ( [69| ) will be minimized and therefore 
the right hand side term in ( [69] ) will be minimized over all the existing server allocation policies. 
In other words, for MW policy and any arbitrary server allocation policy A we have 

D^^(t + l)<iVAL. + (Mirf 



+ 2E 



X{t) A'it + 1)-J2^s [Cs ® /^^(t + 1)) 1^ I X{t) 



< NAl^^ + {MKf 



+ 2E 



X{t) A^{t + 1) - J] TT, {a ® I^{t + 1)) ll I X{t) 



(70) 



where I^^{t + 1) and I^{t + 1) are the allocation matrices of policies MW and A, respectively 
and D^^{t + 1) is the Lyapunov drift of policy MW at time slot t+1. 
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It is important to note that the set of inequahties in (33) determines an open poly tope V 
for which we have V = V — bound{V). In fact, Theorem [2] states that if the vector E[A{t)] 
at each time slot t is strictly inside polytope V then the system is stable. Now, consider an 
MQMS system with arrival processes for which we have E[A{t)] G V at each time slot t. 
Suppose that 5 > is a positive real number such that E[A{t)] + Sljy G hound{V). Therefore, 
E[A{t)\ +5 G P Vt. According to our discussion about randomized policies, there exists a 
randomized policy /^^^ = for which we will have 

T" 



E 



jrnd^^ + 1)) l] 



E[A(t + 1)] + 51 



Therefore, 



E 



x{t) lA'{t + i)-J2^s{Cs®i''{t + i))il] \xit) 



■N 



N 



n=l 



Putting together from ( [69] ) to ( [72| ), we conclude that 

^Mw^^ + 1) < NAl,^ + {MKf - 2S^X^{t) 

n=l 

and according to ([8]) the stability of MQMS system is demonstrated. 



N 



(71) 



-Sj2Xn{t) (72) 



(73) 
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